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LONG  TERM  GOAL 

Our  long  term  goal  is  to  develop  a  tested  model  of  optical  properties  in  the  surf  zone  and 
adjacent  nearshore  ocean,  including  the  influence  of  suspended  sediment,  bubble 
population  and  surface  foam,  to  assess  how  optical  properties  are  related  to  short  tenn 
events  such  as  individual  breaking  wave  crests,  and  to  determine  how  wave-driven  surf 
zone  circulation  influences  spatial  distribution  of  optical  properties  just  offshore  of 
breaking  through  seaward  transport  of  fine  sediment  and  small  persistent  bubbles. 

OBJECTIVES 

The  ability  to  make  optically-based  observations  in  nearshore  waters  is  strongly 
influenced  by  the  presence  of  suspended  sediment  particles  and  of  bubbles,  both  of  which 
are  present  due  to  the  action  of  breaking  waves.  Wave  breaking  is  instrumental  in 
injecting  large  volumes  of  air  into  the  water  column.  This  air  volume  subsequently 
evolves  into  a  distribution  of  bubble  sizes  which  interact  with  the  fluid  turbulence  and  are 
advected  by  the  organized  flow.  Degassing  of  the  water  column  can  additionally  generate 
a  persistent  foam  layer  which  can  prevent  any  optical  penetration  of  the  water  column. 

Our  goal  is  to  develop  time-resolved  models  for  these  processes  in  order  to  make 
predictions  of  optical  properties  of  the  water  column.  To  date,  we  have  begun  this 
process  by  incorporating  a  continuum  description  of  bubble  populations  and  associated 
dynamics  in  2-D  (Ripple)  and  3-D  (Truchas)  Navier-Stokes  solvers.  In  the  continuation 
of  this  work,  our  objectives  are  to: 

1)  Further  develop  the  3-D  Navier  Stokes  framework  to  incorporate  sediment 
particle  phases  as  well  as  chemistry  and  biology  phases  with  possible  interactions 
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with  bubble  populations,  and  to  incorporate  a  foam  generation  model  during 
degassing  processes  at  the  water  surface. 


2)  Further  refine  the  NHWAVE  framework  to  incorporate  similar  extensions,  with  a 
focus  on  surf  zone  scale  (100  m  to  km)  modeling  of  wave  driven  currents  and 
bubble  transport. 

3)  Perfonn  tests  of  both  frameworks  against  recent  field  data  sets. 

APPROACH 

Our  approach  to  the  problem  has  followed  along  two  tracks:  (1)  incorporation  of  a 
comprehensive  model  of  bubble  physics  within  a  3D  LES  hydrodynamic  code.,  for 
application  in  detailed  process  studies  and  development  of  parameterizations  for  use  at 
larger  scale,  and  (2)  development  of  a  model  with  lower  resolution  but  retaining  a  fully  3- 
D  structure,  for  use  in  sur  zone  simulations  at  spatial  scales  on  the  order  of  a  kilometer  or 
so.  The  physics  is  represented  by  a  multiphase  continuum  model,  using  the  fonnalism 
described  by  Drew  and  Passman  (1999),  with  the  details  of  a  fonnulation  for  an  air  and 
water  mixture  found  in  Carrica  et  al  (1999).  In  the  present  project,  we  have  implemented 
a  model  combining  a  water  phase,  a  bubble  phase  with  multiple  bubble  size  (or,  more 
accurately,  mass)  bins.  The  existing  3-D  model  Truchas  has  been  extended  to  include 
Carrica  et  al.’s  polydisperse  multiphase  model  and  an  LES  turbulence  closure  model,  and 
has  been  tested  against  a  number  of  detailed  experimental  data  sets  on  flow  field 
turbulence  and  coherent  vortical  structures  (Ting,  2008;  Ting  and  Nelson,  2011).  The 
effects  of  three-dimensional  obliquely  descending  eddies  (Nadaoka  et  al.,  1989)  and 
downburst  (Kubo  and  Sunamura,  2001)  on  bubble  transport  are  being  investigated  using 
the  3-D  multiphase  model. 

A  framework  for  performing  simulations  of  bubble  fields  and  surface  foam  distributions 
at  kilometer  length  scales  is  being  developed  based  on  a  fully  3-D  nonhydrostatic  wave 
resolving  model.  This  model  will  be  extended  to  incorporate  bubble  phases,  and  will  be 
coupled  to  a  model  for  foam  generation,  transport  and  decay  on  the  water  surface. 

WORK  COMPLETED 

Following  on  work  completed  in  FY11  and  FY12  we  have  implemented  a  dynamic 
turbulence  closure  in  the  full  3-D  VOF  model  for  polydisperse  multiphase  flow  and  have 
carefully  examined  bubble  entrainment  and  liquid-bubble  interaction  under  an  isolated 
breaking  event.  We  have  documented  this  work  in  a  submitted  manuscript  (Derakhti  and 
Kirby,  2013). 

Turning  to  the  problem  of  providing  more  useful  wave-resolved  model  predictions  at 
full  surf  zone  scale,  we  first  developed  a  robust  shock-capturing  scheme  based  on  the 
finite  volume  TVD  approach.  Initially,  we  planned  to  base  the  model  development  on  a 
Boussinesq-type,  depth-integrated  formulation  in  which  the  information  about  void 
fraction  distribution  would  be  incorporated  in  functional  fonn  with  dependence  on  initial 


conditions,  vertical  elevation  and  time  since  the  initiation  of  breaking  event.  This 
approach  let  to  the  development  of  the  Boussinesq  code  FUNWAVE-TVD  (Shi  et  al, 
2012a)  which  provides  a  robust  model  of  a  range  of  surfzone  processes  including 
breaking  and  runup  without  application  of  additional  filtering  steps.  During  this  stage  of 
the  work,  we  decided,  however,  to  concentrate  our  effort  on  a  fully  three-dimensional 
model  which  would  eliminate  the  need  to  impose  the  vertical  structure  of  void  fraction. 
This  led  to  the  development  of  the  code  NHWAVE  (Ma  et  al,  2012a),  which  is  a  shock 
capturing,  fully  nonhydrostatic,  wave  resolving  model  in  surface  and  terrain  following 
sigma  coordinates.  Extension  of  this  model  in  the  context  of  this  project  is  described 
below. 

A  framework  for  performing  simulations  of  bubble  fields  and  surface  foam  distributions 
at  kilometer  length  scales  have  been  developed  based  on  a  fully  3-D  nonhydrostatic  wave 
resolving  model  NHWAVE  (Ma  et  al.  2012a).  The  air  bubble  phase  was  implemented  in 
NHWAVE  following  a  multiphase  description  of  polydisperse  bubble  population  applied 
in  a  3D  VOF  model  by  Ma  et  al.  (2012b).  The  foam  layer  model  is  based  on  a  shallow 
water  fonnulation  with  a  balance  of  drag  forces  due  to  wind  and  water  column  motion 
(Shi  et  al.  2012b).  Foam  mass  conservation  includes  source  and  sink  tenns  representing 
outgassing  of  the  water  column,  direct  foam  generation  due  to  surface  agitation,  and 
erosion  due  to  bubble  bursting.  The  coupled  bubble  and  foam  model  has  been  tested  at 
field-scale  at  FRF,  Duck,  NC. 

RESULTS 

In  addition  to  the  application  of  the  full  3-D  VOF  model  for  surf  zone  breaking  waves 
during  FY1 1  and  FY12,  we  have  carefully  examined  the  model  capabilities  and  accuracy 
to  predict  bubble-entrainment  and  transport  as  well  as  liquid-bubble  interaction  under  an 
isolated  breaking  event.  Figure  1  shows  snapshots  of  the  free  surface  evolution  for  a 
strong  plunging  breaker.  The  model  captures  overturning  jet  impact,  splash-up  process 
and  formation  of  a  bore-like  region.  Finger  shape  structures  can  be  seen  in  the  forward 
splash  (t*  =  0.05  ~  0.3)  and,  comparing  to  Lamarre  and  Melville  (1991,  figure  2a)  the 
time  and  location  t*  —  0.25, x*  =  0.5)  at  which  the  forward  splash  hits  the  undisturbed 
free  surface  is  very  accurately  captured  by  the  model.  During  cavity  fonnation 
(£*  =  —0.15  ~  0)  the  model  predicts  void  fraction  up  to  10%  at  the  jet  toe,  consistent 
with  the  measurement  of  Blenkinsopp  &  Chaplin  (2007,  figure  4,  a  -  c).  The  cavity 
entrapped  by  the  jet  entrains  a  considerable  volume  of  air  during  t*  =  0  ~  0.25.  A 
backward  splash  is  fonned  between  t*  —  0.2  and  t*  —  0.35  and  entrains  some  volume 
of  air.  During  t*  —  0.25  ~  0.5  the  entrained  cavity  collapses  and  big  bubbles  degas  very 
quickly.  The  semicircular  primary  cloud  initially  advances  approximately  with  the  phase 
speed,  but  after  t*  =  0.5  its  horizontal  centroid  becomes  constant  and  then  moves 
backward  slightly  after  t*  =  1.1.  The  secondary  cloud  is  generated  by  the  impact  of  the 
forward  splash  during  t*  —  0.25  to  t*  —  0.65  where  at  t*  —  0.65  we  have  the 
maximum  entrainment  by  a  jet  like  impact  similar  to  the  primary  jet.  The  mean  void 
fraction  of  the  dispersed  bubbles  becomes  more  than  30%  at  t*  —  0.65  and  then 
decreases  gradually  to  1%  at  t*  —  1.7.  While  the  large  bubbles  outgas  very  quickly,  small 
bubbles  are  preferentially  entrained  into  the  coherent  vortices  generated  during  breaking 


and  transported  vertically  by  turbulent  motions,  and  may  remain  in  the  water  column  for 
a  very  long  time.  Bubbles  entrained  by  a  plunging  breaker  can  be  divided  into  three 
different  clouds.  Figure  2  shows  the  3D  dispersed  bubble  plume  evolution,  in  which  the 
two  semicircular  clouds  are  related  to  the  two  downbursts  of  turbulent  motion  under  the 
first  impacting  jet  and  forward  splash,  and  the  third  cloud  represents  bubbles  entrained  by 
the  bore  which  are  transported  by  vortices  behind  the  bore. 

Liquid-bubble  interaction,  i.e.  dispersed  bubble  effects  on  mean  and  turbulent  motions,  is 
still  an  open  question.  In  LES,  the  transport  equation  for  a  resolved  TKE  can  be  obtained 
based  on  the  resolved  velocity  field.  SGS  dissipation  contains  both  shear-  and  bubble- 
induced  dissipation  and  typically  is  much  bigger  than  the  viscous  dissipation.  In  the  case 
of  a  two  phase  flow  with  a  dilute  regime  (a  ~  1),  common  practice  is  to  use  the 
conventional  single  phase  TKE  transport  equation  with  an  additional  term  due  to  a 
correlation  between  fluctuating  concentration  and  vertical  turbulent  velocity  component, 
-pgof  w/  .  High  turbulent  and  dissipative  regions  are  collocated  with  high  void  fraction 
regions.  Due  to  large  void  fraction  (>  10%)  during  active  breaking  period,  the  dilute 
assumption  is  not  valid  anymore,  and  concentration  fluctuations  cannot  be  ignored. 

Figure  3  displays  various  integral  properties  of  the  breakier  for  the  case  of  a  strong 
plunging  breaker.  Figures  3a  and  b  show  the  total  production  rate  by  mean  shear  and 
SGS  dissipation  for  the  simulation  with  and  without  consideration  of  the  dispersed 
bubbles.  Presence  of  the  dispersed  bubbles  reduces  turbulent  mean  shear  production 
while  enhances  the  turbulent  dissipation.  The  total  production  by  buoyancy  is  an  order  of 
magnitude  smaller  than  the  production  by  mean  shear  and  the  dispersed  bubbles  (not 
shown  here).  Figure  3c  shows  the  dispersed  bubbles  damped  TKE  about  20%  to  30%. 
Exception  is  in  0.1  <  t*  <  0.5,  where  TKE  has  not  been  damped  and  even  is  enhanced 
in  the  large  plunging  case.  This  can  be  explained  through  the  work  done  by  the  dispersed 
bubbles  on  the  turbulent  motions.  Figure  3d  shows  total  viscous  (er )  and  SGS 
dissipation  (esgs)  per  unit  crest  width  in  the  breaking  region.  In  all  breaker  types,  most  of 
the  energy  is  dissipated  during  the  first  period  after  breaking,  at  which  bubble-induced 
dissipation  is  about  50%  bigger  than  the  shear-induced  dissipation.  Corresponding 
simulations  without  the  inclusion  of  dispersed  bubbles  predict  smaller  esgs,  about  35%, 
which  is  surprisingly  invariant  with  respect  to  the  different  breaker  types  and  intensity. 
Lamarre  and  Melville  (1991)  approximated  the  total  dissipation  by  estimating  the  total 
energy  flux  difference  between  upstream  and  downstream  of  the  breaking  region.  Table  1 
summarizes  the  dissipation  contributions  predicted  by  the  model  and  the  estimation  from 
Lamarre  and  Melville  (1991).  Comparing  to  the  experimental  estimations,  we  can 
conclude  that  inclusion  of  bubble-induced  dissipation  improves  the  model  results  and  that 
simulations  without  inclusion  of  dispersed  bubbles  underpredict  the  total  dissipation  by 
about  40%  both  in  plunging  and  spilling  breaking. 

The  numerical  codes  FUNWAVE-TVD  (Shi  et  al,  2012a)  and  NHWAVE  (Ma  et  al, 
2012a)  have  been  published  and  are  available  as  open  source  code  for  community  use. 
FUNWAVE-TVD  has  not  been  extended  yet  for  multi-phase  effects,  as  it  was  decided  to 
concentrate  this  effort  on  the  3-D  model  NHWAVE.  Ma  et  al  (2012a)  describe  basic 
testing  of  NHWAVE  as  a  surf  zone  model.  The  model  has  also  been  extended  to 


incorporate  the  multiphase  bubble  model.  At  present,  we  are  performing  more  extensive 
tests  of  the  model’s  ability  to  predict  circulation  and  bubble  void  fraction  over  complex 
topography.  The  NHWAVE  model  is  has  also  been  coupled  to  a  model  for  foam 
distribution  and  transport  on  the  water  surface  (Shi  et  al,  2012b). 

Field-scale  model  testing  is  based  on  an  extensive  history  of  imagery  obtained  at  the 
USAGE  FRF  at  Duck,  NC,  including  ARGUS  imagery  in  the  visible  band  and  IR  band 
imagery  collected  during  the  2010  Surf  Zone  Optics  Experiment  (SZO).  To  incorporate 
AGGUS  measurements  into  the  model  tests,  we  applied  wave  conditions  measured  in 
several  selected  ARGUS  recording  periods.  Figure  4  shows  predictions  of  void  fraction 
(left  panel)  predicted  by  the  NHWAVE  bubble  model  and  foam  coverage  (middle  panel) 
predicted  by  the  foam  layer  model.  IR-based  PIV  surface  velocity  map  shown  in  the  right 
panel  reveals  a  significant  rip  current  between  850-900  m  which  is  also  predicted  by  the 
numerical  model. 

Figure  5  shows  the  time  stack  of  foam  thickness  from  the  model  (left)  compared  with  the 
time  stack  of  nearshore  water  surface  in  infrared  band  from  IR  imagery  (right).  IR 
imagery  is  likely  to  separate  active  breaking  from  passive  foam  in  tenns  of  more  rapidly 
cooling  on  passive  foam.  The  predicted  foam  thickness  looks  consistent  with  foam 
signals  in  IR  imagery,  although  further  work  needs  to  be  done  to  establish  a  quantitative 
correspondence. 


IMPACT/APPLICATIONS 

The  work  proposed  here  would  provide  a  general  framework  for  modeling  the  combined 
effects  of  bubble  distribution,  sediment  load  and  foam  coverage  on  optical  properties  in 
the  surf  zone.  The  model  framework  for  bubble  population  is  intended  to  be  general  in 
nature  and  will  be  applied  at  a  later  date  in  more  computationally  intensive  studies  of 
processes  in  individual  breaking  wave  crests  in  a  wide  range  of  water  depths. 

RELATED  PROJECTS 

Development  of  the  multiphase  aspects  of  the  NHWAVE  model  code  continues  with 
additional  support  from  NSF-PO  through  the  project  “The  interaction  of  waves,  tidal 
currents  and  river  outflows  and  their  effects  on  the  delivery  and  resuspension  of 
sediments  in  the  near  field. 
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Figure  1 :  Snapshots  of  the  time  sequences  of  free  surface  evolution  in  the  breaking  region 
for  a  strong  plunging  breaker. 


Figure  2:  Snapshots  of  3D  bubble  plume  (Isosurface  of  a  —  0.05%)  evolution  in  the 
breaking  region  for  a  strong  plunging  breaker,  (left):  side  view  of  the  3d  results;  (right): 
top  view  of  the  3d  results. 


Figure  3:  Integral  results  for  a  strong  plunging  breaker  case,  a)  The  two  phase  total 
production  rate  by  mean  shear  (solid  line)  simulations  with  the  inclusion  of  the  dispersed 
bubbles  and  (dotted  line)  simulations  without  dispersed  bubbles.  The  reference  value  is 
pL^CcTc1,  where  Lc,  Cc  and  Tc  are  the  wavelength,  phase  speed  and  period  of  the  center 
frequency  of  the  wave  packet,  b)  The  two  phase  total  SGS  dissipation  rate  by  mean  shear 
from  (solid  line)  simulations  with  the  inclusion  of  the  dispersed  bubbles  and  (dotted  line) 
simulations  without  dispersed  bubbles.  The  reference  value  is  pL^CcT^1.  c)  Nonnalized 
total  resolved  TKE  in  the  breaking  region,  (solid  line)  simulations  with  the  inclusion  of 
dispersed  bubbles  and  (dashed  line)  simulations  without  dispersed  bubbles.  The  reference 
value  is  L2cC2S45,  where  S  is  the  wave  packet  steepness,  d)  Total  dissipation  per  unit 
crest  width  in  the  breaking  region  (J/m).  (Thick  solid  line)  Total  dissipation  in  the 
breaking  region,  E^ota[;  (solid  line)  total  dissipation  from  the  simulations  without  the 
inclusion  of  dispersed  bubbles,  E^otai  ■  (dashed-dotted  line)  bubble-induced  dissipation 
EggS  (dashed  line)  shear-induced  dissipation 


x(m) 


Figure  4:  Air  void  fraction  (left)  and  foam  coverage  (middle)  predicted  by  the  coupled 
bubble-foam  model.  IR-based  PIV  surface  velocity  map  (right  panel,  courtesy  of  Andy 
Jessup,  Chris  Chickadel)  reveals  a  significant  rip  current  between  850-900  m  which  is 
also  predicted  by  the  numerical  model. 


Figure  5:  Time  stack  of  foam  thickness  from  the  model  (left)  compared  with  time  stack  of 
nearshore  water  surface  in  infrared  band  from  IR  imagery  (right,  courtesy  of  Andy 
Jessup,  Chris  Chickadel). 


Case  no. 

LM  (J/m) 

etotai  (J/m) 

e total /^tolal  (%) 

tsgs/Ztotal  (%) 

tsqs/ttotal  (%) 

PI 

17.8 

14.7 

63.7 

45.9 

52.9 

P2 

8.6 

7.7 

64.8 

45.4 

53.0 

SP1 

4.3 

2.6 

65.9 

43.2 

53.7 

Table  1:  Total  dissipation  in  the  breaking  region,  Etotai ■  £sgs  and  £sgs  represent  the 
shear-  and  bubble-induced  dissipation,  respectively.  Etotai  is  the  total  dissipation  from  the 
simulations  without  the  inclusion  of  dispersed  bubbles. 


